home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_10_04 / 1004022a < prev    next >
Text File  |  1989-12-31  |  2KB  |  62 lines
  1. /* Listing 1
  2. ** bc program to calculate Chebyshef economized polynomial
  3. ** for evaluation of sin(x) */
  4. /* use bc -l to get c() and s() functions */
  5. define t(x) {                /* sin(x)/x */
  6.     if(x==0)return(1.); /* derivative of s function */
  7.     return (s(x) / x); /* put function to be fit here */ }
  8. define b(x) {
  9.     if (x < 0) return (-x);
  10.     return (x); }
  11. define m(x, y) {
  12.     if (x > y) return (x);
  13.     return (y); }
  14. n = 22; /* number of Chebyshef terms */
  15. scale = 40;
  16. p = a(1.) * 4; /* pi */
  17. b = p * .5; /* upper end of curve fit interval */
  18. a = -b; /* lower end of interval */
  19. /* chebft adapted from Press Flannery et al */
  20. /* "Numerical Recipes" FORTRAN version */
  21. for (k = 1; k <= n; ++k) {
  22. c[k] = 0;
  23. f[k] = t(c((k - .5) * p / n) * (b - a) * .5 + (b + a) * .5);
  24. }
  25. /* because of symmetry, even c[] are 0 */
  26. for (j = 1; j <= n; j += 2) {
  27.     s = 0;
  28.     q = (j - 1) * p / n;
  29.     for (k = 1; k <= n; ++k) s += c(q * (k - .5)) * f[k];
  30.     (c[j] = 2 / n * s); }
  31. /* skip even terms, which are 0 */
  32. for (n = 5; n <= 19; n += 2) {
  33.  /* chebpc */
  34.     for (j = 1; j <= n; ++j) d[j] =  e[j] = 0;
  35.     d[1] = c[n];
  36.     for (j = n - 1; j >= 2; --j) {
  37.     for (k = n - j + 1; k >= 2; --k) {
  38.         s = d[k];
  39.         d[k] = d[k - 1] * 2 - e[k];
  40.         e[k] = s; }
  41.     s = d[1];
  42.     d[1] = c[j] - e[1];
  43.     e[1] = s; }
  44.     for (j = n; j >= 2; --j) d[j] = d[j - 1] - e[j];
  45.     d[1] = c[1] * .5 - e[1];
  46.  /* pcshft */
  47.     g = 2 / (b - a);
  48.     for (j = 2; j <= n; ++j) {
  49.     d[j] *= g;
  50.     g *= 2 / (b - a); }
  51.     for (j = 1; j < n; ++j) {
  52.     h = d[n];
  53.     for (k = n - 1; k >= j; --k) {
  54.         h = d[k] - (a + b) * .5 * h;
  55.         d[k] = h; }
  56.     }
  57.     "Chebyshev Sin fit |x|<Pi/2 coefficients"
  58.     " Maximum Rel Error:"
  59.     m(b(c[n + 2]), b(c[2])) / t(b);
  60.     for (i = 1; i <= n; i += 2)  d[i]; 
  61. }
  62.